;; LTDL_LIBRARY_PATH=/opt/local/lib
;; GUILE_LOAD_PATH=.

;;(import (srfi srfi-1))

(import (rnrs)
        (grid-gen common)
        (grid-gen symath)
)

(module-on real-8)

;;(module-on OpenCL)
;;(module-on OpenCL-CPU)

;;(optimizer:set-max-const-rank 1)

(module-on fortran77)
        
(module-on HD)
        
;;(module-on Lax-Friedrichs)

(module-on Roe)
(module-on Einfeldt)

;;(module-on Harten2)

(module-on Osher)

(module-on MPI 8 8)
;;(module-on OpenMP 2)

;;(module-on OpenMP)

(module-on grid)
(module-on maxima)
(module-on OpenDX)

(module-on multigrid)

(set-gam 1.01)

;; (Harten2:set-eps! (comp-float "eps" 1.0))

(def-float tim 0.0)
(def-float tau)

(def-float cour  0.1)

(define dims '(200 400))
(define dims1 '(400 200))
(define dims0 '(400 400))

(def-float mstar `(* 1.0 ,ASTRO:Msun))
(def-float rstar `(* 1.0 ,ASTRO:Rsun))
(def-float rdisk `(* 5.0 ,ASTRO:Rsun))


(def-float xw `(* 20.0 ,ASTRO:Rsun))
(def-float yw `(* 20.0 ,ASTRO:Rsun))

(def-float xmin1 `(* -0.5 ,xw))
(def-float xmax1 `(- 0.0 (* 0.25 (/ ,xw ,(- (car dims) 1.0)))))
(def-float xmin2 `(+ 0.0 (* 0.25 (/ ,xw ,(- (car dims) 1.0)))))
(def-float xmax2 `(*  0.5 ,xw))

(def-float ymin `(* -0.5 ,yw))
(def-float ymax `(*  0.5 ,yw))

(def-float Pmax `(/ (* 2.0 ,MATH:Pi) ,(ASTRO:OmegaK mstar rdisk)))


(def-float maxtemp (CGS 1.0e4 'K))
(def-float maxv `(* 10.0 ,rstar ,(ASTRO:OmegaK mstar rstar)))

(def-float tmax `(* ,Pmax 1.0))
(def-float rho_disk (CGS 1.0e-13 'rho))
(def-float t_disk   (CGS 1.0e4 'K))

(set-scal 'cm rdisk)
(set-scal 's `(* ,Pmax 1.0e-5))
;;(set-scal 'g ASTRO:Msun)
(set-scal 'rho (2CGS rho_disk 'rho))

(def-float rhofon `(* ,rho_disk 1.0e-5))
(def-float tfon   t_disk)

(define p_disk `(* ,rho_disk ,PHYS:Rgas ,t_disk))
(def-float pfon `(* ,rhofon ,PHYS:Rgas ,tfon))

(define grdf0 (lambda (x y z)
  (list `(+ ,xmin1 (* (- ,xmax2 ,xmin1) ,x))
        `(+ ,ymin (* (- ,ymax ,ymin) ,y))
        z
  )
))

(define grdf1 (lambda (x y z)
  (list `(+ ,xmin1 (* (- ,xmax1 ,xmin1) ,x))
        `(+ ,ymin (* (- ,ymax ,ymin) ,y))
        z
  )
))

(define grdf2 (lambda (x y z)
  (list `(+ ,xmin2 (* (- ,xmax2 ,xmin2) ,x))
        `(+ ,ymin (* (- ,ymax ,ymin) ,y))
        z
  )
))

(define grdf-circ (lambda (rmin rmax gi ng dms)
  (lambda (r0 ph0 z)
    (let* ((ny (cadr dms))
           (r     `(+ ,rmin (* ,r0 (- ,rmax ,rmin))))
           (dphi  `(/ (* 2 ,MATH:Pi) ,ng))
           (dphi1 `(/ ,dphi ,ny))
           (phi0  `(+ (* ,dphi ,gi) (* 0.5 ,dphi1)))
           (dphi2 `(- ,dphi ,dphi1))
           (phi  `(+ ,phi0 (* ,ph0 ,dphi2)))
          )
      (list `(* ,r (cos ,phi))
            `(* ,r (sin ,phi))
             z
      )
    )
  )
))


;(define mgr (multigrid ((gr1 (grid:cart dims grdf1))
;                        (gr2 (grid:cart dims grdf2))
;                       )
;              (grid:join gr1 'x+ '(y) gr2 'x- '(y))
;            )
;)

(def-float rmax `(* ,rdisk 1.5))
(def-float rmin `(* ,rstar 0.9))

;(define mgr (multigrid ((gr1 (grid:cart dims1 (grdf-circ rmin rmax 0 2 dims1)))
;                        (gr2 (grid:cart dims1 (grdf-circ rmin rmax 1 2 dims1)))
;                       )
;              (grid:join gr1 'y+ '(x) gr2 'y- '(x))
;              (grid:join gr2 'y+ '(x) gr1 'y- '(x))
;            )
;)

(define mgr (multigrid ((gr1 (grid:cart dims0 (grdf-circ rmin rmax 0 1 dims0))))
              (grid:join gr1 'y+ '(x) gr1 'y- '(x))
            )
)

;(define mgr (grid:cart dims0 grdf0))


(define star-cond (lambda (gr i)
  `(<= (+ (** ,(gr 'x i) 2) (** ,(gr 'y i) 2)) (** ,rstar 2.0))
))

(mgr 'add-bnd 'free-in star-cond
  (lambda (gr i)
    (beg
      (set-val (gr 'r i) rhofon)
      (set-val (gr 'p i) pfon)
      (set-val (gr 'u i) 0.0)
      (set-val (gr 'v i) 0.0)
    )
  )
)


;;;;(mgr 'add-bnd `x- 'free-in)

(define kepl (lambda (x y)
  (let* ((r  `(sqrt (+ (** ,x 2) (** ,y 2))))
         (v  (ASTRO:OmegaK mstar r))
         (vx `(* -1.0 ,y ,v))
         (vy `(* ,x ,v))
        )
    `(,vx ,vy 0.0)
;;    `((* -1.0 ,vx) (* -1.0 ,vy) 0.0)
  )
))


;(gr 'add-corr-pred 'tmax (CGS 1.0e6 'K))
;(gr 'add-corr-pred 'rhomax `(* 1.0e3 ,rho_disk))
;(gr 'add-corr-pred 'vmax `(* ,(ASTRO:VK mstar rstar) 10.0))

;(gr 'add-corr 'smooth)

;(gr 'add-corr 'old)

(mgr 'add-corr
    (lambda (gr i)
      (beg
        (set-val (gr 'r1 i) rhofon)
        (set-val (gr 'p1 i) pfon)
        (set-val (gr 'u1 i) 0.0)
        (set-val (gr 'v1 i) 0.0)
      )
    )
)

(mgr 'add-corr 'err)

(set! PHYS:Potential (lambda (x y z)
  `(/ (* -1.0 ,PHYS:G ,mstar) (sqrt (+ (** ,x 2) (** ,y 2))))
))


(define tP (lambda (t) `(/ ,t ,Pmax)))

(define nrm (lambda (x) `(/ ,x ,ASTRO:Rsun)))

(def-float dvc)

(define chop-temp (lambda (gr)
  (lambda (vr)
    (beg
    ; `(if (> ,(vr 'temp) maxtemp)
    ;     ,(set-val (gr 'p) `(* ,(gr 'r) ,PHYS:Rgas ,maxtemp))
    ;     ,nop
    ;  )
      (set-val (gr 'p) `(* ,(gr 'r) ,PHYS:Rgas ,t_disk))
     `(if (> ,(vr 'u2) (** ,maxv 2.0))
         ,(beg
            (set-val dvc `(sqrt (/ (** ,maxv 2.0) ,(vr 'u2))))
            (set-val (gr 'u) `(* ,dvc ,(gr 'u)))
            (set-val (gr 'v) `(* ,dvc ,(gr 'v)))
          )
         ,nop
      )
    )
  )
))

(define save-dx
  (beg
    (OpenDX:save mgr "disk.dx" 
     `(("positions" (x y) ,(2CGS-lam 'cm))
       ("rho" r ,(2CGS-lam '(/ g (** cm 3.0))))
       ("p" p ,(2CGS-lam '(/ dyn (** cm 2.0))))
       ("v" (u v) ,(2CGS-lam '(/ cm s)))
      )
    )
  )
)

(def-float nt 0.0)

(def-float dt)
(def-float dtm)
(def-float ta 0.0)
(def-float t0)
(def-float t1)

(define rate `(floor (* (/ (- ,tmax ,tim) (/ ,tim ,nt)) (/ ,ta ,nt))))

(define init-disk (lambda (gr i)
  `(if (and (<= (+ (** ,(gr 'x i) 2) (** ,(gr 'y i) 2)) (** ,rdisk 2.0))
            (>= (+ (** ,(gr 'x i) 2) (** ,(gr 'y i) 2)) (** ,rstar 2.0))
       )
      ,(beg
         (set-val (gr 'r i) rho_disk)
         (set-val (gr 'p i) p_disk)
         (set-val (gr 'u i) (car  (kepl (gr 'x i) (gr 'y i))))
         (set-val (gr 'v i) (cadr (kepl (gr 'x i) (gr 'y i))))
       )
      ,(beg
         (set-val (gr 'r i) rhofon)
         (set-val (gr 'p i) pfon)
         (set-val (gr 'u i) 0.0)
         (set-val (gr 'v i) 0.0)
       )
   )
))

(define save-vr
  (beg
    (to-file "grid.dat" 'main 'replace
      (apply beg
        (map (lambda (gr)
               (beg
                 (gr 'forall-grid-d '(0 1) '(x y)
                       (lambda (x y)
                         (fprint (cur-ptr 0) (cur-ptr 1) x y)
                       )
                 )
                 (fprint "-----------")
               )
             )
             (mgr 'subgrids)
        )
      )
    )
  )
)

(code-gen "test"
  (beg
;;    (timer-tick)
    (set-val t0 (time))
    (mgr 'init-grid)
    
    save-vr
    
    (mgr 'init init-disk)
;;;    (do-while (and (< tim tmax) (= (gr 'err) 0))
    (do-while (< tim tmax)
      (timer-tick)
      (mgr 'update-bnd)
      
      (do-every 30 'sec
        save-dx
      )

      (if (= (mgr 'err) 0)
        (beg
          (set-val (mgr 'tau) (- tmax tim))
          (mgr 'CFL cour)
        )
      )

      (measure-time dt
        (mgr 'load-grid-vals1)
        (mgr 'advection)
        (mgr 'correction)
        (mgr 'update-grid-vals)
      )

;;      (mgr 'forall-grid-local '(var) (chop-temp mgr))
      
      (if (> nt 0)
          (set-val dtm (min dt dtm))
          (set-val dtm dt)
      )
      
      (set-val ta (+ ta dt))
      (set-val nt (+ nt 1.0))

      (do-every 1 'sec
        (set-val t1 (time))
        (main-p
          (print "t=" (tP tim)
                 "tau=" (tP (mgr 'tau))
                 (floor (* (/ (- tmax tim) tim) (- t1 t0)))
          )
        )
      )
 
      (if (= (mgr 'err) 0)
          (set-val tim (+ tim (mgr 'tau)))
          (beg
            (set-val (mgr 'tau) (* 0.5 (mgr 'tau)))
            ;;(main-p (do-every 1 'sec (print "new Tau=" (tP (mgr 'tau)))))
          )
      )
    )
    (mgr 'update-bnd)
    save-dx
  )
)

(display "\n")
